home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / astronomy / skyview / !SkyView / c / Planets < prev    next >
Encoding:
Text File  |  1993-08-26  |  32.5 KB  |  927 lines
  1. /********************************************************/
  2. /*              Planets Module for SkyView              */
  3. /*                                                      */
  4. /*                  (c)1992 N P Hawkes                  */
  5. /*                                                      */
  6. /*                   Displays Planets                   */
  7. /********************************************************/
  8.  
  9. #include "menu.h"
  10. #include "dbox.h"
  11. #include "bbc.h"
  12. #include "wimpt.h"
  13. #include "dbox.h"
  14. #include "sprite.h"
  15. #include "res.h"
  16. #include "resspr.h"
  17. #include "werr.h"
  18.  
  19. #include "sv_header.h"
  20. #include "planets.h"
  21. #include "datime.h"
  22. #include "radec.h"
  23. #include "ecliptic.h"
  24. #include "riset.h"
  25.  
  26. #include <stdio.h>
  27. #include <stdlib.h>
  28. #include <string.h>
  29. #include <math.h>
  30.  
  31. /********************************************************/
  32. /*                     Constants                        */
  33. /********************************************************/
  34. #define DISP_NAME "Planets"   /* Entry in Display menu. */
  35. #define SEL_NAME  "Planet"    /* Entry in Select menu.  */
  36.  
  37. #define MAX_PLANS 10          /* Max no. of planets.    */
  38. #define FILE_NAME "PlanData"  /* Name of data file.     */
  39.  
  40. #define NAME_LEN  15          /* Length of name field.  */
  41. #define NCOEFFS    4          /* # of coeffs in polynoms*/
  42.  
  43. #define SPR_NAME "planet"     /* Name of sprite.        */
  44. #define SPR_X0 -12            /* Bounding box of        */
  45. #define SPR_Y0 -12            /* sprite, OS units,      */
  46. #define SPR_X1  14            /* relative to            */
  47. #define SPR_Y1  14            /* plotting position.     */
  48. #define SPR_COL  5            /* Sprite x offset.       */
  49. #define SPR_ROW  4            /* Sprite y offset.       */
  50. #define SPR_MODE MODE_20      /* Sprite mode.           */
  51.  
  52. /*  Set effective altitude of horizon for calculations  */
  53. /* of rising and setting.                               */
  54. #define HORALT (REAL)0.0088
  55.  
  56. /* Convergence criteria for times of phenomena:         */
  57. #define MAX_ITER 5
  58. #define CONVERGE 5
  59.  
  60. /********************************************************/
  61. /*                New types of variable.                */
  62. /********************************************************/
  63. /* Structure describing one planet.                     */
  64. typedef struct {
  65.   char name[NAME_LEN+1];  /* Name of planet.            */
  66.   REAL meanl[NCOEFFS];    /* Coeffs for mean longitude. */
  67.   REAL lperi[NCOEFFS];    /*   " for longit. of perihel.*/
  68.   REAL eccen[NCOEFFS];    /*   " for eccentricity.      */
  69.   REAL inclin[NCOEFFS];   /*   " for inclination.       */
  70.   REAL lanode[NCOEFFS];   /*   " for long. of asc. node.*/
  71.   REAL semimaj;           /* Semi-major axis.           */
  72.   REAL angdiam;           /* Angular diameter at 1 AU.  */
  73.   REAL stanmag;           /* Standard magnitude.        */
  74.   REAL sdist;             /* Apparent dist from Sun.    */
  75.   REAL edist;             /* Apparent dist from Earth.  */
  76.   REAL phase;             /* Phase (1 = full).          */
  77.   REAL ra;                /* Right Ascension of planet. */
  78.   REAL dec;               /* Declination of planet.     */
  79.   int  horiz_id;          /* id in Horiz window.        */
  80.   int  vert_id;           /* id in Vert window.         */
  81. } planstr;
  82.  
  83. /* Structure containing ecliptic rectangular            */
  84. /* coordinates and velocity components.                 */
  85. typedef struct {
  86.   double x;
  87.   double y;
  88.   double z;
  89.   double xdot;
  90.   double ydot;
  91.   double zdot;
  92. } rect_coords;
  93.  
  94. /********************************************************/
  95. /*                  Global Variables                    */
  96. /********************************************************/
  97. static BOOL display_flag = TRUE;   /* 'Enabled' flag.   */
  98. static BOOL data_valid   = FALSE;  /* Validity of list. */
  99. static int moduleid;               /* Module ID.        */
  100.  
  101. static FILE *fileptr;              /* Ptr to data file. */
  102. static char fmt[256];              /* Format string.    */
  103. static planstr plan_data[MAX_PLANS];/* Main data store. */
  104. static int plan_count = 0;         /* No. of planets.   */
  105.  
  106. static menu plan_menu;    /* For selecting a planet.    */
  107. static sprite_id spr_id;  /* ID of planet sprite.       */
  108.  
  109. static double T_now;      /* Centuries from 0/1/1900.   */
  110.  
  111. /********************************************************/
  112. /*                 Function Prototypes                  */
  113. /********************************************************/
  114. static void plans_buildfn(void);
  115. static void plans_selectfn(selectfn_reasoncode reason);
  116. static void selection_details(void);
  117. static BOOL cul_calc(int id, REAL altmin,
  118.                  REAL *altptr, REAL *azimptr, int *hourptr, int *minptr);
  119. static BOOL rs_calc(int id, BOOL whichcalc,
  120.                  REAL *altptr, REAL *azimptr, int *hourptr, int *minptr);
  121. static BOOL plans_displayfn(BOOL *enabptr);
  122. static os_error *plans_plotfn(int x, int y, int id);
  123. static void plans_infofn(int id);
  124. static BOOL read_plandata(planstr *pptr);
  125. static char *trim(char *str);
  126. static void plan_altaz_etc(int id, rect_coords *earthptr,
  127.                            REAL *altptr, REAL *azimptr);
  128. static void earth_coords(double T, rect_coords *earthptr);
  129. static void plan_coords(int id, double T, rect_coords *planptr);
  130. static double ang_element(REAL coeff[NCOEFFS], double T);
  131. static double ang_elemen1(REAL coeff[NCOEFFS], double T);
  132. static double element(REAL coeff[NCOEFFS], double T);
  133. static void plan_geo_elatlong(rect_coords *planptr, rect_coords *earthptr,
  134.                               REAL *elatitptr, REAL *elongitptr);
  135. static void plan_geo_elatlong_etc(
  136.                               rect_coords *planptr, rect_coords *earthptr, 
  137.                               REAL *elatitptr, REAL *elongitptr,
  138.                               REAL *sdistptr,  REAL *edistptr,
  139.                               REAL *phaseptr);
  140. static void plan_radec(int id, double jul, REAL *raptr, REAL *decptr);
  141. static void plan_radecfn(int id, observerstr *ob_ptr, int hour, int min,
  142.                          REAL *raptr, REAL *decptr);
  143. static void plan_altaz_any(int id, observerstr *ob_ptr, int hour, int min,
  144.                            REAL *altptr, REAL *azimptr);
  145.  
  146. /********************************************************/
  147. /*                Initialisation Function               */
  148. /********************************************************/
  149. BOOL plans_initfn(int moduleno, modulestr *plans)
  150. {
  151.   int i;
  152.   sprite_id nametype_id;
  153.   os_error *errptr;
  154.  
  155. /* Record the module number of this module.             */
  156.   moduleid = moduleno;
  157.  
  158. /* Open data file.                                      */
  159.   fileptr = res_openfile(FILE_NAME, "r");
  160.   if (fileptr==NULL) return FALSE;
  161.  
  162. /* Build format string for reading name of planet.      */
  163.   sprintf(fmt, " %%%ic", NAME_LEN);
  164.  
  165. /* Read data, one planet at a time.  Stop on first      */
  166. /* error (probably EOF).  Stay within array bounds.     */
  167.   while (plan_count < MAX_PLANS && \
  168.          read_plandata(&plan_data[plan_count]))
  169.     plan_count++;
  170.  
  171. /* Close file.  Fatal error if it won't close.          */
  172.   if (fclose(fileptr) != 0)
  173.     werr(FATAL, "Planets Error: Can't close Data file");
  174.  
  175. /* Quit if no items have been read.                     */
  176.   if (plan_count == 0) return FALSE;
  177.  
  178. /* Build menu for Select function.                      */
  179. /* Use menu_new for first item, menu_extend thereafter. */
  180.   plan_menu = menu_new(SEL_NAME ":", plan_data[0].name);
  181.   if (plan_menu == NULL) return FALSE;
  182.   for (i=1; i<plan_count; i++)
  183.     menu_extend(plan_menu, plan_data[i].name);
  184.  
  185. /* Fill in the fields of the modulestr.                 */
  186.   plans->buildfn  = plans_buildfn;
  187.   plans->selectfn = plans_selectfn;
  188.   plans->dispfn   = plans_displayfn;
  189.   plans->infofn   = plans_infofn;
  190.   plans->initial  = display_flag;
  191.   plans->display_entry = DISP_NAME;
  192.   plans->display_menu  = NULL;
  193.   plans->select_entry  = SEL_NAME;
  194.   plans->select_menu   = plan_menu;
  195.  
  196. /* Set up the sprite_id structure which identifies the  */
  197. /* planet sprite.                                       */
  198. /* First set up a name-type identifier.                 */
  199.   nametype_id.s.name = SPR_NAME;
  200.   nametype_id.tag    = sprite_id_name;
  201. /* Then get the address of the sprite, and put it in    */
  202. /* the global sprite identifier spr_id.                 */
  203.   errptr = sprite_select_rp(resspr_area(), &nametype_id, &spr_id.s.addr);
  204.   if (errptr != NULL) return FALSE;
  205. /* Finally tag the identifier spr_id as address-type.  */
  206.   spr_id.tag = sprite_id_addr;
  207.  
  208.   return TRUE;
  209. }
  210.  
  211. /*------------------------------------------------------*/
  212. /*     Function to read in the data for one planet.     */
  213. /* Data consist of:                                     */
  214. /*   Planet's name,                                     */
  215. /*   Coeffs for polynomials for five orbital elements,  */
  216. /*   Semi-major axis of orbit,                          */
  217. /*   Angular diameter at 1 AU,                          */
  218. /*   Standard magnitude.                                */
  219. /* All must be present.  Function returns TRUE if all   */
  220. /* data read OK, FALSE otherwise.                       */
  221. /* A suitable source of data is Duffett-Smith p133.     */
  222. /*------------------------------------------------------*/
  223. static BOOL read_plandata(planstr *pptr)
  224. {
  225. /* Temp storage for polynomial coefficients:            */
  226.   float c0, c1, c2, c3;
  227. /* Temp storage for other floating-point values:        */
  228.   float f0, f1, f2;
  229.  
  230. /* Read name.                                           */
  231.   if ( fscanf(fileptr, fmt, pptr->name) != 1 )
  232.     return FALSE;
  233.  
  234. /* Put a terminating \0 into the name string.  Then     */
  235. /* trim trailing blanks.                                */
  236.   pptr->name[NAME_LEN] = '\0';
  237.   trim(pptr->name);
  238.  
  239. /* Read coeffs for mean longitude.                      */
  240.   if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
  241.     return FALSE;
  242.   pptr->meanl[0] = CONV * (REAL)c0;
  243.   pptr->meanl[1] =        (REAL)c1;
  244.   pptr->meanl[2] = CONV * (REAL)c2;
  245.   pptr->meanl[3] = CONV * (REAL)c3;
  246.  
  247. /* Read coeffs for longitude of perihelion.             */
  248.   if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
  249.     return FALSE;
  250.   pptr->lperi[0] = CONV * (REAL)c0;
  251.   pptr->lperi[1] = CONV * (REAL)c1;
  252.   pptr->lperi[2] = CONV * (REAL)c2;
  253.   pptr->lperi[3] = CONV * (REAL)c3;
  254.  
  255. /* Read coeffs for eccentricity.                        */
  256.   if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
  257.     return FALSE;
  258.   pptr->eccen[0] = (REAL)c0;
  259.   pptr->eccen[1] = (REAL)c1;
  260.   pptr->eccen[2] = (REAL)c2;
  261.   pptr->eccen[3] = (REAL)c3;
  262.  
  263. /* Read coeffs for inclination.                         */
  264.   if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
  265.     return FALSE;
  266.   pptr->inclin[0] = CONV * (REAL)c0;
  267.   pptr->inclin[1] = CONV * (REAL)c1;
  268.   pptr->inclin[2] = CONV * (REAL)c2;
  269.   pptr->inclin[3] = CONV * (REAL)c3;
  270.  
  271. /* Read coeffs for longitude of ascending node.         */
  272.   if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
  273.     return FALSE;
  274.   pptr->lanode[0] = CONV * (REAL)c0;
  275.   pptr->lanode[1] = CONV * (REAL)c1;
  276.   pptr->lanode[2] = CONV * (REAL)c2;
  277.   pptr->lanode[3] = CONV * (REAL)c3;
  278.  
  279. /* Read semimajor axis, angular diam & standard magn.   */
  280.   if ( fscanf(fileptr, "%f %f %f", &f0, &f1, &f2) != 3 )
  281.     return FALSE;
  282.   pptr->semimaj = (REAL)f0;
  283.   pptr->angdiam = (REAL)f1;
  284.   pptr->stanmag = (REAL)f2;
  285.  
  286.   return TRUE;
  287. }
  288.  
  289. /*------------------------------------------------------*/
  290. /*  Function which trims trailing blanks from a string. */
  291. /*------------------------------------------------------*/
  292. static char *trim(char *str)
  293. {
  294.   char *ptr;
  295.  
  296.   for (ptr = str+strlen(str); ptr>str && *(ptr-1)==' '; ptr--)
  297.     ;
  298.   *ptr = '\0';
  299.  
  300.   return str;
  301. }
  302. /********************************************************/
  303. /*                 List-Building Function               */
  304. /********************************************************/
  305. static void plans_buildfn(void)
  306. {
  307.   int i;                /* Planet ID.                   */
  308.   plotobj planet;       /* Plotobj under construction.  */
  309.   rect_coords earth;    /* Current coords of Earth.     */
  310.  
  311.   /* Bounding box of plotting symbol, relative to       */
  312.   /* position of symbol:                                */
  313.   static wimp_box size = {SPR_X0, SPR_Y0, SPR_X1, SPR_Y1};
  314.  
  315. /* If module is disabled, do not rebuild list.  Instead,*/
  316. /* set validity flag to FALSE and return.               */
  317.   if (!display_flag)
  318.   {
  319.     data_valid = FALSE;
  320.     return;
  321.   }
  322.  
  323. /* Calculate current heliocentric ecliptic rectangular  */
  324. /* coords of Earth.                                     */
  325.   T_now = ob_data.jul/36525.0;
  326.   earth_coords(T_now, &earth);
  327.  
  328. /* For each planet in turn:                             */
  329.   for (i=0; i<plan_count; i++)
  330.   {
  331.     /* Build plotobj, and record interesting planetary  */
  332.     /* data for possible use by info function.          */
  333.     planet.id = i;
  334.     planet.module = moduleid;
  335.     planet.plotfn  = plans_plotfn;
  336.     plan_altaz_etc(i, &earth, 
  337.                    &planet.alt, &planet.azim);
  338.     planet.size = size;
  339.  
  340.     /* Offer this to the main program.                  */
  341.     addobj(planet, &plan_data[i].horiz_id, &plan_data[i].vert_id);
  342.   }
  343.  
  344. /* Set flag to show that data in list is now valid.     */
  345.   data_valid = TRUE;
  346.  
  347.   return;
  348. }
  349.  
  350. /*------------------------------------------------------*/
  351. /*  Function to calculate current alt & azim (and other */
  352. /*  data of interest) for planet No. 'id'.  T_now and   */
  353. /*coords of Earth (pointed to by earthptr) must be valid*/
  354. /*------------------------------------------------------*/
  355. static void plan_altaz_etc(int id, rect_coords *earthptr,
  356.                            REAL *altptr, REAL *azimptr)
  357. {
  358.   rect_coords planet;
  359.   REAL elatit, elongit;
  360.   REAL ra, dec;
  361.  
  362. /*Get planet's heliocentric rectangular ecliptic coords.*/
  363.   plan_coords(id, T_now, &planet);
  364.  
  365. /* Derive geocentric ecliptic latit. & longit, plus     */
  366. /* other planetary data of interest.                    */
  367.   plan_geo_elatlong_etc(&planet, earthptr,
  368.            &elatit, &elongit,
  369.            &plan_data[id].sdist, &plan_data[id].edist, &plan_data[id].phase);
  370.  
  371. /* Convert ecliptic lat & long to RA & Dec.             */
  372.   ecliptic_radec(elatit, elongit, T_now, &ra, &dec);
  373.  
  374. /* Convert RA and Dec to altitude & azimuth.            */
  375.   radec_altaz(ra, dec, &ob_data, ob_data.sid, altptr, azimptr);
  376.  
  377.   plan_data[id].ra  = ra;
  378.   plan_data[id].dec = dec;
  379.  
  380.   return;
  381. }
  382.  
  383. /*------------------------------------------------------*/
  384. /* Function to calculate the RA and Dec of planet 'id'  */
  385. /* at the instant implied by 'jul'.                     */
  386. /*------------------------------------------------------*/
  387. static void plan_radec(int id, double jul, REAL *raptr, REAL *decptr)
  388. {
  389.   double T;
  390.   REAL elatit, elongit;
  391.   rect_coords earth, planet;
  392.  
  393.   T = jul/36525.0;
  394.  
  395. /* Get Earth's heliocentric rectangular ecliptic coords.*/
  396.   earth_coords(T, &earth);
  397.  
  398. /* Get same for planet.                                 */
  399.   plan_coords(id, T, &planet);
  400.  
  401. /* Derive geocentric ecliptic latit. and longit.        */
  402.   plan_geo_elatlong(&planet, &earth, &elatit, &elongit);
  403.  
  404. /* Convert latit. and longit. to RA and Dec.            */
  405.   ecliptic_radec(elatit, elongit, T, raptr, decptr);
  406.  
  407.   return;
  408. }
  409.  
  410. /*------------------------------------------------------*/
  411. /* Function to calculate the RA and Dec of planet 'id'  */
  412. /* at a specified instant.  Same as plan_radec above,   */
  413. /* but arguments conform to requirements of 'Riset' fns.*/
  414. /*------------------------------------------------------*/
  415. static void plan_radecfn(int id, observerstr *ob_ptr, int hour, int min,
  416.                          REAL *raptr, REAL *decptr)
  417. {
  418.   plan_radec(id, datime_julian(ob_ptr, hour, min), raptr, decptr);
  419.  
  420.   return;
  421. }
  422.  
  423. /*------------------------------------------------------*/
  424. /* Function to calculate the alt and azim of planet 'id'*/
  425. /* at any time on the day designated in the specified   */
  426. /* observerstr (which also defines the location).       */
  427. /*------------------------------------------------------*/
  428. static void plan_altaz_any(int id, observerstr *ob_ptr, int hour, int min,
  429.                              REAL *altptr, REAL *azimptr)
  430. {
  431.   double jul;      /* Days since noon GMT 0/1/1900.     */
  432.   REAL   sid;      /* Corresponding local siderial time.*/
  433.   REAL ra, dec;
  434.  
  435.   jul = datime_julian(ob_ptr, hour, min);
  436.   sid = datime_siderial(ob_ptr, jul, hour, min);
  437.  
  438. /* Get RA and Dec of planet.                            */
  439.   plan_radec(id, jul, &ra, &dec);
  440.  
  441. /* Convert to altitude and azimuth.                     */
  442.   radec_altaz(ra, dec, ob_ptr, sid, altptr, azimptr);
  443.  
  444.   return;
  445. }
  446.  
  447. /*------------------------------------------------------*/
  448. /*  Function to find rectangular heliocentric ecliptic  */
  449. /*     coords and velocity components for the Earth.    */
  450. /* Assumes semimajor axis is 1.0 AU, and inclination is */
  451. /*  zero.  Method based on that in NAO Technical Note   */
  452. /*      No.46, Section 5 (Dec 1978 / Feb 1981), p5.     */
  453. /*------------------------------------------------------*/
  454. static void earth_coords(double T, rect_coords *ptr)
  455. {
  456.   double L;      /* Mean longitude.                     */
  457.   double M;      /* Mean anomaly.                       */
  458.   double ec;     /* Eccentricity.                       */
  459.   double E;      /* Eccentric anomaly.                  */
  460.   double WW;     /* Longitude of perihelion.            */
  461.   double x, y, r, u, v;
  462.   double K = 0.0172021;
  463.  
  464. /* Get elements of Sun's "orbit".                       */
  465.   ecliptic_sunorbit(T, &L, &M, &ec);
  466.  
  467. /* Convert mean longit. of Sun to that of Earth.        */
  468.   L += DblPI;
  469.  
  470. /* Get longitude of perihelion.                         */
  471.   WW = L - M;
  472.  
  473. /* Solve Kepler's eqn for eccentric anomaly.            */
  474.   ecliptic_kepler(M, ec, &E);
  475.  
  476. /* Get coords and velocity components.                  */
  477.   x = cos(E) - ec;
  478.   y = sqrt(1.0 - ec*ec) * sin(E);
  479.   r = sqrt(x*x + y*y);
  480.   u = -(K * sin(E))/r;
  481.   v =  (K * sqrt(1.0 - ec*ec) * cos(E))/r;
  482.  
  483.   ptr->x    = x*cos(WW) - y*sin(WW);
  484.   ptr->y    = x*sin(WW) + y*cos(WW);
  485.   ptr->z    = 0.0;
  486.   ptr->xdot = u*cos(WW) - v*sin(WW);
  487.   ptr->ydot = u*sin(WW) + v*cos(WW);
  488.   ptr->zdot = 0.0;
  489.  
  490.   return;
  491. }
  492. /*------------------------------------------------------*/
  493. /*  Function to find rectangular heliocentric ecliptic  */
  494. /*    coords and velocity components for planet 'id'.   */
  495. /*  Method based on that in NAO Technical Note No.46,   */
  496. /*        Section 5 (Dec 1978 / Feb 1981), p5.          */
  497. /*------------------------------------------------------*/
  498. static void plan_coords(int id, double T, rect_coords *ptr)
  499. {
  500.   double L;      /* Mean longitude.                     */
  501.   double WW;     /* Longitude of perihelion.            */
  502.   double ec;     /* Eccentricity.                       */
  503.   double i;      /* Inclination.                        */
  504.   double O;      /* Longitude of ascending node.        */
  505.   double a;      /* Semi major axis.                    */
  506.   double M;      /* Mean anomaly.                       */
  507.   double w;      /* Argument of the perihelion.         */
  508.   double E;      /* Eccentric anomaly.                  */
  509.   double x,  y,  r,  u,  v;
  510.   double x1, y1,     u1, v1;
  511.   double cosw, sinw;
  512.   double cosi, sini;
  513.   double cosO, sinO;
  514.   double K = 0.0172021;
  515.  
  516. /* Get elements of orbit at time T.                     */
  517.   L = ang_elemen1(plan_data[id].meanl, T);
  518.   WW= ang_element(plan_data[id].lperi, T);
  519.   ec=     element(plan_data[id].eccen, T);
  520.   i = ang_element(plan_data[id].inclin,T);
  521.   O = ang_element(plan_data[id].lanode,T);
  522.   a = (double)plan_data[id].semimaj;
  523.  
  524.   M = L - WW;
  525.   w = WW - O;
  526.  
  527. /* Solve Kepler's eqn for eccentric anomaly.            */
  528.   ecliptic_kepler(M, ec, &E);
  529.  
  530. /* Get coords and velocity components.                  */
  531.   x = a * (cos(E) - ec);
  532.   y = a * sqrt(1.0 - ec*ec) * sin(E);
  533.   r = sqrt(x*x + y*y);
  534.   u = -(K * sqrt(a) * sin(E))/r;
  535.   v =  (K * sqrt(a*(1.0 - ec*ec)) * cos(E))/r;
  536.  
  537.   cosw = cos(w);  cosi = cos(i);  cosO = cos(O);
  538.   sinw = sin(w);  sini = sin(i);  sinO = sin(O);
  539.  
  540.   x1 = x*cosw - y*sinw;
  541.   y1 = x*sinw + y*cosw;
  542.   u1 = u*cosw - v*sinw;
  543.   v1 = u*sinw + v*cosw;
  544.  
  545.   ptr->x    = x1*cosO - y1*sinO*cosi;
  546.   ptr->y    = x1*sinO + y1*cosO*cosi;
  547.   ptr->z    = y1*sini;
  548.   ptr->xdot = u1*cosO - v1*sinO*cosi;
  549.   ptr->ydot = u1*sinO + v1*cosO*cosi;
  550.   ptr->zdot = v1*sini;
  551.  
  552.   return;
  553. }
  554.  
  555. /*------------------------------------------------------*/
  556. /*   Function to calculate an angular orbital element   */
  557. /*   from its polynomial coefficients.  Result is in    */
  558. /*   radians and in the range - 2 pi < element < 2 pi.  */
  559. /*------------------------------------------------------*/
  560. static double ang_element(REAL coeff[NCOEFFS], double T)
  561. {
  562.   return fmod(element(coeff, T), DblTWO_PI);
  563. }
  564.  
  565. /*------------------------------------------------------*/
  566. /*   Function to calculate an orbital element from its  */
  567. /*              polynomial coefficients.                */
  568. /*------------------------------------------------------*/
  569. static double element(REAL coeff[NCOEFFS], double T)
  570. {
  571.   double el;
  572.   double T2 = T * T;
  573.  
  574.   el  = (double)coeff[0];
  575.   el += (double)coeff[1] * T;
  576.   el += (double)coeff[2] * T2;
  577.   el += (double)coeff[3] * T2 * T;
  578.  
  579.   return el;
  580. }
  581.  
  582. /*------------------------------------------------------*/
  583. /* Function to calculate an orbital element from coeffs */
  584. /*      intended for use with a modified T term.        */
  585. /* See Duffett-Smith p.132.                             */
  586. /*------------------------------------------------------*/
  587. static double ang_elemen1(REAL coeff[NCOEFFS], double T)
  588. {
  589.   double B;
  590.   double el;
  591.   double T2 = T * T;
  592.  
  593. /* Subtract integer rotations from T term, for better   */
  594. /* accuracy.                                            */
  595.   el  = (double)coeff[1] * T;
  596.   B   = DblTWO_PI*(el - floor(el));
  597.  
  598.   el  = (double)coeff[0];
  599.   el += B;
  600.   el += (double)coeff[2] * T2;
  601.   el += (double)coeff[3] * T2 * T;
  602.  
  603.   return fmod(el, DblTWO_PI);
  604. }
  605.  
  606. /*------------------------------------------------------*/
  607. /* Function to calculate geocentric ecliptic lat & long */
  608. /* for a planet, plus other planetary data of interest. */
  609. /*            Allows for light travel time.             */
  610. /*------------------------------------------------------*/
  611. static void plan_geo_elatlong_etc(
  612.                            rect_coords *planptr,
  613.                            rect_coords *earthptr,
  614.                            REAL *elatptr,  REAL *elongptr,
  615.                            REAL *sdistptr, REAL *edistptr, REAL *phaseptr)
  616. {
  617.   double R;           /* True distance from Earth.      */
  618.   double rho;         /* Apparent distance from Earth.  */
  619.   double xp, yp, zp;  /* Apparent geo. rectang. coords. */
  620.   double tau;         /* Light travel time.             */
  621.   double elongit;     /* App. geocent. ecliptic long.   */
  622.   double helongit;    /* App. heliocent. ecliptic long. */
  623.   double phase;       /* App. phase of planet.          */
  624.   double x, y, z;
  625.   double xdot, ydot, zdot;
  626.  
  627. /* Get true distance from Earth.                        */
  628.   x = planptr->x - earthptr->x;
  629.   y = planptr->y - earthptr->y;
  630.   z = planptr->z - earthptr->z;
  631.   R = sqrt(x*x + y*y + z*z);
  632.  
  633. /* Get light travel time.                               */
  634.   tau = R / 173.14;
  635.  
  636. /* Get apparent distance from Earth.                    */
  637.   xdot = planptr->xdot - earthptr->xdot;
  638.   ydot = planptr->ydot - earthptr->ydot;
  639.   zdot = planptr->zdot - earthptr->zdot;
  640.   xp = x - tau*xdot;
  641.   yp = y - tau*ydot;
  642.   zp = z - tau*zdot;
  643.   rho = sqrt(xp*xp + yp*yp + zp*zp);
  644.   *edistptr = (REAL)rho;
  645.  
  646. /* Get apparent geocentric ecliptic lat & long.         */
  647.    elongit  = atan2(yp, xp);
  648.   *elatptr  = (REAL)asin(zp / rho);
  649.   *elongptr = (REAL)elongit;
  650.  
  651. /* Get apparent distance from Sun.                      */
  652.   xp = planptr->x  -  tau * planptr->xdot;
  653.   yp = planptr->y  -  tau * planptr->ydot;
  654.   zp = planptr->z  -  tau * planptr->zdot;
  655.   *sdistptr = (REAL)sqrt(xp*xp + yp*yp + zp*zp);
  656.  
  657. /* Get apparent heliocentric ecliptic longitude.        */
  658.   helongit = atan2(yp, xp);
  659.  
  660. /* Get phase (see Duffett-Smith p.139).                 */
  661.   phase = 0.5 * (1.0 + cos(elongit - helongit));
  662.   if (phase < 0.0) phase = 0.0;
  663.   if (phase > 1.0) phase = 1.0;
  664.   *phaseptr = (REAL)phase;
  665.  
  666.   return;
  667. }
  668.  
  669. /*------------------------------------------------------*/
  670. /* Function to calculate geocentric ecliptic lat & long */
  671. /* for a planet, allowing for light travel time.        */
  672. /*------------------------------------------------------*/
  673. static void plan_geo_elatlong(
  674.                            rect_coords *planptr,
  675.                            rect_coords *earthptr,
  676.                            REAL *elatptr,  REAL *elongptr)
  677. {
  678.   double R;           /* True distance from Earth.      */
  679.   double rho;         /* Apparent distance from Earth.  */
  680.   double xp, yp, zp;  /* Apparent geo. rectang. coords. */
  681.   double tau;         /* Light travel time.             */
  682.   double elongit;     /* App. geocent. ecliptic long.   */
  683.   double x, y, z;
  684.   double xdot, ydot, zdot;
  685.  
  686. /* Get true distance from Earth.                        */
  687.   x = planptr->x - earthptr->x;
  688.   y = planptr->y - earthptr->y;
  689.   z = planptr->z - earthptr->z;
  690.   R = sqrt(x*x + y*y + z*z);
  691.  
  692. /* Get light travel time.                               */
  693.   tau = R / 173.14;
  694.  
  695. /* Get apparent distance from Earth.                    */
  696.   xdot = planptr->xdot - earthptr->xdot;
  697.   ydot = planptr->ydot - earthptr->ydot;
  698.   zdot = planptr->zdot - earthptr->zdot;
  699.   xp = x - tau*xdot;
  700.   yp = y - tau*ydot;
  701.   zp = z - tau*zdot;
  702.   rho = sqrt(xp*xp + yp*yp + zp*zp);
  703.  
  704. /* Get apparent geocentric ecliptic lat & long.         */
  705.    elongit  = atan2(yp, xp);
  706.   *elatptr  = (REAL)asin(zp / rho);
  707.   *elongptr = (REAL)elongit;
  708.  
  709.   return;
  710. }
  711.  
  712.  
  713. /********************************************************/
  714. /*              Object-Selecting Function               */
  715. /********************************************************/
  716. static void plans_selectfn(selectfn_reasoncode reason)
  717. {
  718.   int menu_hit;
  719.  
  720.   switch (reason)
  721.   {
  722.     case new_selection:
  723.       /* Look at list of menu hits to find which        */
  724.       /* planet is to be selected.                      */
  725.       menu_hit = module_submenu_hits[0];
  726.  
  727.       /* If no hits on Planet menu, flag the failure    */
  728.       /* and return.  Also return if hit is off the end.*/
  729.       if (menu_hit <= 0 || menu_hit > plan_count)
  730.       {
  731.         selection.selected_OK = FALSE;
  732.         return;
  733.       }
  734.  
  735.       selection.selected_OK = TRUE;
  736.       /* Fill in Selection details for requested        */
  737.       /* planet.  Allow for menu entries being          */
  738.       /* indexed from 1 instead of 0.                   */
  739.       selection.id = menu_hit - 1;
  740.       selection_details();
  741.       break;
  742.  
  743.     case window_selection:
  744.       /* A selection has been made by clicking on a     */
  745.       /* planet symbol in a window.                     */
  746.     case recalculate_data:
  747.       /* Observer details have changed, and plotting    */
  748.       /* lists have been rebuilt.                       */
  749.       /* */
  750.       /*In either case, the planet ID is in selection.id*/
  751.       selection_details();
  752.       break;
  753.  
  754.     case timeonly_recalculate:
  755.       /* As recalculate_data, but only the time of day   */
  756.       /* (and possibly the direction of view - this is   */
  757.       /* of no interest) have changed.                   */
  758.       selection.horiz_id    = plan_data[selection.id].horiz_id;
  759.       selection.vert_id     = plan_data[selection.id].vert_id;
  760.       /* Decide if planet can be seen now.               */
  761.       selection.now      = selection.horiz_id != NOWHERE || \
  762.                            selection.vert_id  != NOWHERE;
  763.       break;
  764.  
  765.     case sel_info_request:
  766.       /* Write info into text buffer.                   */
  767.       /* Use the Info function.                         */
  768.       plans_infofn(selection.id);
  769.       break;
  770.  
  771.     default:
  772.       /* Unknown option.  Do nothing.                   */
  773.       break;
  774.  
  775.     }
  776.  
  777.   return;
  778. }
  779.  
  780.  
  781. /*------------------------------------------------------*/
  782. /*  Function to fill in the details of rising, setting  */
  783. /*                   and culminating.                   */
  784. /*------------------------------------------------------*/
  785. static void selection_details(void)
  786. {
  787.   REAL rise_alt;  /*Altitude at calculated rising time. */
  788.   REAL set_alt;   /*Altitude at calculated setting time.*/
  789.  
  790.   selection.horiz_id = plan_data[selection.id].horiz_id;
  791.   selection.vert_id  = plan_data[selection.id].vert_id;
  792.   selection.now      = selection.horiz_id != NOWHERE || \
  793.                        selection.vert_id  != NOWHERE;
  794.  
  795.   selection.culminating = cul_calc(selection.id, HORALT,
  796.                               &selection.cul_alt,
  797.                               &selection.cul_azim,
  798.                               &selection.cul_hour,
  799.                               &selection.cul_min);
  800.   selection.rising      = rs_calc(selection.id, RISECALC,
  801.                               &rise_alt,
  802.                               &selection.rise_azim,
  803.                               &selection.rise_hour,
  804.                               &selection.rise_min);
  805.   selection.setting     = rs_calc(selection.id, SETCALC,
  806.                               &set_alt,
  807.                               &selection.set_azim,
  808.                               &selection.set_hour,
  809.                               &selection.set_min);
  810.   return;
  811. }
  812.  
  813. /*------------------------------------------------------*/
  814. /*      Function to derive details of culmination.      */
  815. /* If planet culminates, fn finds alt, azim, hour & min */
  816. /* and returns TRUE.  Returns FALSE if no culmination.  */
  817. /*------------------------------------------------------*/
  818. static BOOL cul_calc(int id, REAL altmin,
  819.                 REAL *altptr, REAL *azimptr, int *hourptr, int *minptr)
  820. {
  821.   BOOL ok;
  822.  
  823. /* Guess a time of day.  Take midday.                   */
  824.   *hourptr = 12;
  825.   *minptr  =  0;
  826.  
  827. /* Try to improve guessed time.                         */
  828.   ok = riset_cul(
  829.        &ob_data, MAX_ITER, CONVERGE, plan_radecfn, id, hourptr, minptr);
  830.  
  831. /* Give up if time did not converge.                    */
  832.   if (!ok) return FALSE;
  833.  
  834. /* Check that planet would be higher than 'altmin'.     */
  835.   plan_altaz_any(id, &ob_data, *hourptr, *minptr, altptr, azimptr);
  836.   return (*altptr >= altmin);
  837. }
  838.  
  839. /*------------------------------------------------------*/
  840. /*   Function to derive details of rising or setting.   */
  841. /* If phenomenon occurs, fn finds alt, azim, hour & min */
  842. /*      and returns TRUE.  Returns FALSE otherwise.     */
  843. /*------------------------------------------------------*/
  844. static BOOL rs_calc(int id, BOOL whichcalc,
  845.                     REAL *altptr, REAL *azimptr, int *hourptr, int *minptr)
  846. {
  847.   BOOL ok;
  848.  
  849. /* Guess a time of day.  Take midday.                   */
  850.   *hourptr = 12;
  851.   *minptr  =  0;
  852.  
  853. /* Try to improve this guess.                           */
  854.   ok = riset_riset(&ob_data, MAX_ITER, CONVERGE, plan_radecfn, id,
  855.                    whichcalc, HORALT, hourptr, minptr);
  856.  
  857. /* Give up if time did not converge to a valid value.   */
  858.   if (!ok) return FALSE;
  859.  
  860. /* Check that planet is indeed visible at this time.       */
  861.   plan_altaz_any(id, &ob_data, *hourptr, *minptr, altptr, azimptr);
  862.   return (*altptr > ZERO);
  863. }
  864.  
  865.  
  866. /********************************************************/
  867. /*                   Display Function                   */
  868. /********************************************************/
  869. static BOOL plans_displayfn(BOOL *enabptr)
  870. {
  871. /* Toggle Enable/Disable flag.                          */
  872.   display_flag = !display_flag;
  873.  
  874. /* Inform main program of new status.                   */
  875.   *enabptr = display_flag;
  876.  
  877. /* If module is now enabled, but plotting list is not   */
  878. /* valid, call list-building function.                  */
  879.   if (display_flag && !data_valid)
  880.     plans_buildfn();
  881.  
  882. /* Windows will always need updating.                   */
  883.   return TRUE;
  884. }
  885.  
  886.  
  887. /********************************************************/
  888. /*                   Plotting Function                  */
  889. /********************************************************/
  890. static os_error *plans_plotfn(int x, int y, int id)
  891. {
  892.   return sv_plotsprite(&spr_id, SPR_MODE, x, y, SPR_COL, SPR_ROW);
  893. }
  894.  
  895.  
  896. /********************************************************/
  897. /*                     Info Function                    */
  898. /********************************************************/
  899. static void plans_infofn(int id)
  900. {
  901.   char p_text[RADEC_TEXTLEN]; /*Text buff for RA & Dec.*/
  902.   REAL app_diam;  /* Apparent angular diameter.        */
  903.   REAL logarg;
  904.   REAL app_mag;   /* Apparent magnitude (see           */
  905.                   /* Duffett-Smith p.139).             */
  906.  
  907. /* Get RA and Dec in text form.                        */
  908.   radec_text(plan_data[id].ra, plan_data[id].dec, p_text);
  909.  
  910. /* Get apparent angular diameter and magnitude.        */
  911.   app_diam = plan_data[id].angdiam / plan_data[id].edist;
  912.   logarg = (plan_data[id].edist * plan_data[id].sdist) /
  913.            (REAL)sqrt((double)plan_data[id].phase);
  914.   app_mag  = (REAL)5.0*(REAL)log10((double)logarg) + plan_data[id].stanmag;
  915.  
  916.   sprintf(infoptr,
  917.   "%s%s\n%s%.1f\n%s%.1f%s\n%s%.2f%s%.1f%c\n%s",
  918.     "Planet: ", plan_data[id].name,
  919.     "Approx visual mag. ", (double)app_mag,
  920.     "Distance ", (double)plan_data[id].edist, " AU",
  921.     "Phase ", (double)plan_data[id].phase,
  922.     ", Ang diam ", (double)app_diam, '"',
  923.     p_text);
  924.  
  925.   return;
  926. }
  927.